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In primordial gas, molecular hydrogen forms primarily through associative 



O 

detachment of H and H, thereby destroying the H . The H anion can also be 
destroyed by a number of other reactions, most notably by mutual neutraliza- 



tion with protons. However, neither the associative detachment nor the mutual 
neutralization rate coefficients are well determined: both may be uncertain by 
as much as an order of magnitude. This introduces a corresponding uncertainty 
into the H 2 formation rate, which may have cosmological implications. Here, we 
examine the effect that these uncertainties have on the formation of H2 and the 
cooling of protogalactic gas in a variety of situations. We show that the effect is 
particularly large for protogalaxies forming in previously ionized regions, affect- 
ing our predictions of whether or not a given protogalaxy can cool and condense 
within a Hubble time, and altering the strength of the ultraviolet background 
that is required to prevent collapse. 



Subject headings: atomic data — galaxies: formation — molecular data — molec- 
ular processes — stars: formation 
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1. 



Introduction 



It has long been argued that molecular hydrogen, H2, must play a central role in the 
cooling of primordial gas in the first protogalaxies. This is simply because it is the only 
coolant present in significant quantities that remains effective at temperatures below I0 4 K 
(Saslaw & Zipoy 1967; Peebles & Dicke f968; Matsuda, Sato, & Takeda 1969). Subsequent 
numerical work has only served to confirm this (e.g., Hutchins 1976; Palla, Salpeter, & Stahler 
1983; Tegmark et al. 1997; Machacek, Bryan, & Abel 2001; Abel, Bryan, & Norman 2002). 
It is now clear that one of the keys to understanding the earliest episodes of star formation 
in the cosmos is a detailed understanding of the chemistry and thermodynamics of H 2 . 
The chemistry of primordial gas has recently been reviewed by a number of authors (Abel 
et al. 1997; Galli & Palla 1998; Standi, Lepp, & Dalgarno 1998; Lepp, Stancil, & Dalgarno 
2002). These authors find that considerable chemical complexity is possible despite the 
limited number of elements available (essentially only hydrogen, helium, and lithium, plus 
isotopic counterparts such as deuterium or 3 He) . But much of this complexity arises from the 
chemistry of minor coolants such as HD or LiH, or trace molecules such as Hjj" or HeH + . These 
molecules do not play a significant role in the chemistry of H 2 (Abel et al. 1997; Glover 2001) 
and the cooling that they provide is generally unimportant compared to that coming from 
H 2 . It is thus usually unnecessary to model their chemistry (although HD cooling can become 
important in gas with an atomic hydrogen number density Wh > 10 4 cm~ 3 and temperature 
T < 200 K; see Flower & Pineau des Forets 2001 or Nakamura & Umemura 2002). The 
omission of these molecules allows substantial simplifications to be made. In particular, Abel 
et al. (1997) showed that the formation and destruction of H 2 can be accurately followed 
over a wide range of temperatures and densities with as few as 28 reactions. In many 
circumstances, this number can be reduced further. For instance, four of the reactions deal 
with the ionization and recombination of helium, and so play no role in gas in which helium 
remains neutral. Furthermore, three of the photochemical reactions - the photoionization of 
H and H 2 and the photodissociation of H 2 by absorption into the continuum of the Lyman 
and Werner band systems (Allison & Dalgarno 1969) - require photons with energies above 
the Lyman limit and so are important only if a strong source of hard ultraviolet photons is 
present. 

The key reactions responsible for the formation of H 2 are easily summarized. Direct 
radiative association of atomic hydrogen is strongly forbidden (Gould & Salpeter 1963), 
and so the main gas-phase pathway by which H 2 is formed makes use of the H~ ion as an 
intermediary. H~ is formed by the slow radiative association reaction 



H + c 



H- + 7, 




and is then destroyed by a fast associative detachment reaction with atomic hydrogen, form- 
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ing H 2 : 

H- + H^H 2 + e" (2) 
Molecular hydrogen can also be formed by a similar chain of reactions involving : 

H + H+ - H+ + 7 , (3) 
H + H+ -> H 2 + H + , (4) 

but in most circumstances the H~ pathway dominates (Glover 2003) as we now explain. 

Various reactions compete with associative detachment to destroy H~, of which the 
most important are mutual neutralization with protons 

fT + H+^H + H, (5) 

and photodetachment by infrared and/or optical photons 

fr+7^H + e-. (6) 

At very high redshifts (z > 100), the cosmic microwave background (CMB) temperature is 
large enough to produce a substantial photodetachment rate (Galli & Palla 1998; Stancil et 
al. 1998) and so only a few of the H~ ions survive for long enough to form H 2 . At these 
redshifts, H 2 formation is dominated by the slower pathway. However, at lower redshifts, 
photodetachment of H~ by the CMB rapidly becomes unimportant; and at the redshifts 
of interest in this paper its effect is negligible. The net rate of H 2 formation is therefore 
determined by two main factors - the fractional ionization of the gas, which controls the 
rate of the initiating radiative association (reaction 1), and the fraction of the resulting H~ 
that is destroyed by associative detachment (reaction 2). A similar state of affairs exists for 
formation of H 2 by the pathway, with the main alternative destruction mechanism in 
this case being dissociative recombination: 

H+ + e~ -> H + H. (7) 

Once formed, H 2 can be destroyed by charge transfer with H + , 

H 2 + H+ -> H+ + H, (8) 

by collisional dissociation by free electrons or atomic hydrogen 1 

H 2 + e- -> H + H + e", (9) 
H 2 + H -> H + H + H, (10) 



dissociation due to collisions with molecular hydrogen is also possible, but collisions with atomic 
hydrogen dominate due to the fact that the molecular fraction of primordial gas is small (typically, 
n H2 < 5 x 10~ 3 n H ; sec Susa et al. 1998, Nishi & Susa 1999, Oh & Haiman 2002 for a detailed discus- 
sion of why this is the case). 
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or by photodissociation via the Solomon process (Stecher & Williams 1967) 

H 2 + 7 ^H + H. (11) 

In addition to the reactions listed above, several others are required to complete our 
chemical network. Collisional ionization of H by electrons 

H + c~ -> H + + e"+e- (12) 

and radiative recombination of H + 

H+ + e- -> H + 7, (13) 

must be included if we wish to model the evolution of the fractional ionization of the gas 

correctly. The other reactions include the dissociative recombination of (reaction 7 above) 
together with 

Fr + e~ -> H + e^+e", (14) 

FT + H -> H + H + e", (15) 

H" + H + -> H+ + e", (16) 

H+ + 7 -> H + H+. (17) 

These reactions all play a role in regulating either the H~ or the abundance, and so 
affect H 2 formation. However, it should be noted that in most circumstances, the influence 
of these reactions on the H 2 formation rate is relatively small. 

There is considerable variation in the accuracy with which the rates and rate coefficients 
of the reactions in this network are known. Some, such as the radiative association reaction 
that forms H~ (reaction 1), have rate coefficients which have been determined to a high level 
of accuracy (see the discussion in Abel et al. 1997). Others, however, such as the H 2 charge 
transfer reaction (reaction 8) discussed by Savin et al. (2004), have rate coefficients which 
are considerably more uncertain. In this paper, we are concerned with two reactions for 
which a large degree of uncertainty exists: the destruction of H~ by associative detachment 
with H (reaction 2) and by mutual neutralization with H + (reaction 5). In situations where 
the photodetachment of H~ (reaction 6) is unimportant, it is the competition between these 
two reactions that determines what fraction of the H~ formed by reaction 1 goes on to form 
n 2 . So substantial variations in the rate coefficients of these reactions can result in major 
variations in the H 2 formation rate. This is turn may have significant consequences for the 
cooling of primordial gas in situations of cosmological importance. It is therefore important 
to determine the extent to which the existing uncertainties in the associative detachment 
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and mutual neutralization rate coefficients affect the conclusions that we can draw regarding 
the formation of H 2 in primordial gas. 

To investigate this, we have performed several sets of simulations of protogalactic col- 
lapse in which we have varied the two rate coefficients over a range of plausible values, in 
an effort to see how sensitive is the outcome for various initial conditions to the particular 
collision data chosen. In Section 2 of this paper, we review the current state of knowledge re- 
garding the values of the associative detachment and mutual neutralization rate coefficients. 
In Section 3 we describe the code used to perform our simulations. The initial conditions 
used in our simulations are discussed in Section 4, and the results of the simulations are 
presented in Section 5. We conclude with a brief discussion in Section 6. 

2. Data for Mutual Neutralization and Associative Detachment 

Here we review the published theoretical and experimental work for the associative 
detachment reaction 2 and the mutual neutralization reaction 5 at temperatures T and 
collision energies E relevant for cosmology. In specific, we focus on work at T < 10 4 K or 
E < 1 eV. The corresponding published rate coefficients for reactions 2 and 5 are shown in 
Figs. 1 and 2, respectively. 

2.1. Associative Detachment 

The only measurement at cosmological temperatures of the associative detachment re- 
action 2 has been carried out by Schmeltekopf, Fehsenfeld, & Ferguson (1967). They used 
a fast-flowing afterglow system and claimed their measured rate coefficient is accurate to 
within a factor of 2. This was quickly followed with a semi-classical theoretical calculation 
by Dalgarno & Browne (1967) using a simple extension on the theory for radiative associ- 
ation. Browne & Dalgarno (1969) subsequently revised their earlier calculations using new 
potentials for the intermediate H2 molecular anion. They presented results for their poten- 
tials (1+3) and (2+4). Both of these results from Browne & Dalgarno are shown in Fig. 1. 
All of the above theory takes into account the lowest two electronic states of H^", namely the 
attractive X 2 £+ state and the repulsive B 2 E+ state. Bieniek & Dalgarno (1979) revisited 
this reaction ten years later treating the problem as one of scattering by a complex potential 
(i.e., resonant scattering theory) but appear to have only accounted for the X 2 E+ state. All 
of these calculations are in reasonable agreement with the laboratory results, to within the 
quoted experimental factor of 2 uncertainty. 
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Fig. I. — Published experimental and theoretical rate coefficients for the associative de- 
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Fig. 2. — Published theoretical and experimental rate coefficients for the mutual neutraliza- 
tion reaction H~ + H + — > H + H at temperatures relevant for early Universe chemistry. 
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Sakimoto (1989) calculated the rate coefficient for reaction 2, following the quantum 
methodology of Bieniek & Dalgarno (1979). However, Sakimoto paid attention to lower 
collision energies than considered by Bieniek & Dalgarno and also used new results for the 
X 2 £+ potential. The resulting rate coefficient is a factor of almost 4 times larger than the 
measurement of Schmeltekopf et al. (1967). Sakimoto writes that this may be due to his 
having not taken into account the contribution of the B 2 £+ state. Sakimoto also presents 
a classical Langevin limit rate coefficient which falls outside of the experimental uncertainty 
limits. 

Using resonant scattering theory and a local complex potential, Launay, Le Dourneuf, 
& Zeippen (1991) calculated the rate coefficient for reaction 2 taking into account only the 
X 2 £+ state. Launay et al. present results for their potentials VI and V2. Each resulting 
rate coefficient lies within the factor of 2 experimental error bars. However, the energy 
dependence calculated by Launay et al. differs from that of Browne & Dalgarno (1969) 
above ~ 2, 000 K. Also, it appears that the results using the potential V2 are incorrect as is 
discussed in more detail by Cizek, Horacek, & Domcke (1998). 

The most recent calculation for reaction 2 is that of Cizek et al. (1998). They use 
a nonlocal resonance model and a new potential for the X 2 £+ state. They do not take 
into account the B 2 £+ state. The resulting rate coefficient lies a factor of ~ 3 above 
the experimental results of Schmeltekopf et al. (1967). Their classical Langevin limit rate 
coefficient is in good agreement with the experiment but lies a factor of R3 2 below the 
Langevin results of Sakimoto (1989). 

Taking the lower limit for the experimental rate coefficient and the highest calculated 
value, yields almost an order of magnitude spread. This is surprisingly large for such a 
simple reaction. A partial explanation for the disagreement between the various theoretical 
calculations may lie in the Hg potentials used. Amaya-Tapia, Cisneros, & Russek (1986) 
showed that theory had not yet converged to a single set of potentials for either the X 2 £+ 
or B 2 £+ state of H2 . The most recent calculations for reaction 2 use differing potentials 
for the X 2 £+ state and have ignored the B 2 S+ state (Sakimoto 1989; Launay et al. 1991; 
Cizek et al. 1998) and all yield differing results. With theory unable to converge to the same 
value for the rate coefficient, it appears that the only hope of improving our understanding 
of reaction 2 lies with carrying out new laboratory measurements. 
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2.2. Mutual Neutralization 

The first published theoretical rate coefficient for reaction 5 at cosmological tempera- 
tures appears to be the Landau-Zener (LZ) results of Bates & Lewis (1955). Moseley, Aberth, 
& Peterson (1970) are the only group to have measured reaction 5 at collision energies below 
3 eV. They carried out measurements to energies as low as 0.1 eV. Using their cross section 
results and extrapolating below 0.1 eV, they produced an experimentally-derived rate coef- 
ficient (Moseley et al. 1970; Peterson et al. 1971). Their experimental rate coefficient is a 
factor of pa 3 times larger than the LZ results of Bates & Lewis (1955). 

Dalgarno & McCray (1973) write that the rate coefficient for 

IT+X+^H + X (18) 

is "fairly insensitive to the nature of the ions X + " . They recommend a rate coefficient for 
all such reactions which is based on the results of Peterson et al. (1971) for reaction 5. 
Surprisingly, the rate coefficient of Dalgarno & McCray is a factor of ~ 2 times smaller than 
that of Peterson et al. This may be a partial explanation as to why Dalgarno & McCray 
estimate that their recommended rate coefficient is accurate to "perhaps a factor of 2". 
Pradad & Huntress (1980), in their commonly cited compilation of rate coefficients for gas 
phase chemistry, use the recommended rate coefficient of Dalgarno & McCray (1973). This 
then made its way into the recommended data of Duley & Williams (1984) and then into 
Shapiro & Kang (1987). 

Fussen & Kubach (1986) calculated reaction 5 using a quantum close-coupling treat- 
ment. Their cross section results are in good agreement with the LZ calculations of Bates 
& Lewis (1955) but a factor of 3 times smaller than the measurements of Moseley et al. 
(1970). Croft, Dickinson, & Gadea (1999) used the cross section calculations of Fussen & 
Kubach to produce a rate coefficient. 

Dalgarno & Lepp (1987) present two rate coefficients for reaction 5, a high value and a 
low value. The higher value is supposedly based on the results of Moseley et al. (1970), but 
as can be seen in Fig 2, lies a factor of ~ 4 below their results. The lower value is supposed 
to be derived from the experimental results of Szucs et al. (1984) and Peart, Bennett, & 
Dolder (1985) which lie close to the theoretical cross section predictions of Bates & Lewis 
(1955). However, these experimental results were carried out for collision energies in excess 
of 5 eV which is of little relevance for early Universe chemistry. Also, Fig. 2 shows that the 
lower recommended value of Dalgarno & Lepp is a factor of ^ 3 times smaller than that 
of Bates & Lewis (1955). The cause for the differences for each rate coefficient of Dalgarno 
& Lepp and their stated sources is unclear. Perhaps both represent typographical errors. 
The lower rate coefficient of Dalgarno & Lepp was adopted by Abel et al. (1997) for their 
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primordial chemistry. 

Clearly, there is still a large uncertainty in the rate coefficient for reaction 5. Is is true 
that LZ and quantum close-coupling calculations are in good agreement for T < 10 4 K. But 
the only measurement at the relevant collision energies lies a factor of 3 higher. Cross 
section measurements at collision energies above 3 eV have been carried out by a number 
of groups (Szucs et al. 1984; Peart et al. 1985; Peart & Hayton 1992) and these are in 
reasonable agreement with the calculations of Bates & Lewis (1955) and Fussen & Kubach 
(1986). However, since the work of Moseley et al. (1970) no further measurements have 
been carried out at collision energies relevant for early Universe chemistry. Additionally, no 
explanation has been given to date as to why the data of Moseley et al. (1970) may be wrong 
(cf., Szucs et al. 1984; Peart et al. 1985; Peart & Hayton 1992). It is clear that in order 
to resolve this issue, new laboratory measurements are needed for reaction 5 at collision 
energies below 1 eV. 



3. Computational method 

To aid us in assessing the impact of the existing uncertainties in the associative detach- 
ment and mutual neutralization rate coefficients on the cooling and collapse of primordial 
gas, we have performed a number of simulations of the collapse of gas into protogalactic halos 
using a numerical technique known as smoothed particle hydrodynamics, or SPH (Gingold & 
Monaghan 1977; Lucy 1977; Monaghan 1992). SPH is a Lagrangian method for simulating 
astrophysical fluid flows, in which the fluid is represented by an ensemble of 'particles', with 
flow quantities at a particular point obtained by averaging over an appropriate subset of 
neighbouring SPH particles. It should be noted that the particles used in SPH simulations 
in no sense correspond to actual gas particles; in typical astrophysical simulations, the for- 
mer have masses greatly in excess of a solar mass and so represent more than 10 57 actual 
gas particles. Instead, each SPH particle simply corresponds to a small region of the flow, 
with fixed mass but variable volume (cf., Eulerian grid-based codes, such as that described 
in Stone & Norman 1992, where the volume of each grid cell is fixed, but the mass of fluid 
contained within it varies). To perform our simulations, we used a a modified version of the 
Gadget SPH code (Springel, Yoshida, & White 2001). Full details of our modifications are 
given elsewhere (Jappsen et al. 2005b), but we briefly summarize them here. 

The most significant modification that we have made to the code is the addition of the 
necessary framework for following the non-equilibrium chemistry of H2 in the protogalactic 
gas. To incorporate this into the code, we associate a set of chemical abundances with each 
SPH particle. Just as with the other fluid properties, such as density or internal energy, 
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these abundances represent averages over the local fluid flow. For each SPH particle, we 
therefore must solve a set of chemical equations of the form 

^ = C l -D l x l (19) 

where Xi is the fractional abundance of species i, computed with respect to the total number 
density of hydrogen nuclei (i.e., x^ = rii/n, where rii is the number density of species % and n 
is the number density of hydrogen nuclei), and where Cj and Di are terms representing the 
chemical creation and destruction of species %. The values of the creation and destruction 
terms generally depend on both the density and temperature of the gas, as well as on its 
chemical composition. 

In our code, we follow the chemistry of six species: e - , H + , H, H~, H^, and H 2 . Two of 
these species, H~ and H^, reach chemical equilibrium on very short timescales. In the case 
of H~, one can readily show that equilibrium is reached on a timescale 

*eq,H- - (^ad^H + ^1%+)^ , (20) 

where »h is the hydrogen atom density, %+ is the H + density, and k^d and /c mn are the 
associative detachment and mutual neutralization rate coefficients, respectively. We assume 
that other processes responsible for destroying H~, such as photodetachment (reaction 6) can 
be neglected. If %+ <C (& a dAmn)^H, then the associative detachment term dominates and 
since A; ad ~ 10~ 9 cm 3 s -1 (to within an order of magnitude), it follows that £ eq ,H- ~ lO 9 ^ 1 s. 
On the other hand, if n H + (&adAmn)^H, then the mutual neutralization term dominates, 
in which case t cqiH - ~5x 10 5 T^ 2 n R l s, again to within an order of magnitude. In either 
case, the equilibrium timescale is very short. The largest value we find for it occurs within 
the neutral, low density gas present at the beginning of runs A1-A9 (see Section 5.1 below), 
but even here, where riu ~ 10~ 3 , its value is only t e q,H- ~ 10 12 s, which is significantly smaller 
than either the cooling time and the free-fall time of the gas (which are both typically much 
longer than a Myr, as we demonstrate in Section 5.3). In more ionized gas, or at higher 
densities, £ e q,H- becomes even smaller relative to the cooling time t coo \ and the free-fall time 
t s . 

The case of is very similar: one can show that it reaches chemical equilibrium on a 
timescale 

*eq,H+ - ( k ctnH + k^ll^y 1 , (21) 

where k ct and k& T are the rate coefficients of the charge transfer and dissociative recombi- 
nation reactions which are the main processes responsible for removing H^" from the gas. 
(For reference, these are reactions 4 and 7 in Table 1). Evaluating this expression, we obtain 
t cqH + ~ 1.6 x lO 9 ^ 1 s in the low ionization limit, or t eqH + ~ lO 8 ^^ s in the high ionization 
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limit. Again, these timescales are both very much shorter than either the cooling time or 
the free-fall time of the gas. 

Because H~ and reach equilibrium so quickly, we do not attempt to follow their 
chemical equilibrium directly in our simulations. Instead, their abundances are computed 
only as required, under the assumption that both reach equilibrium instantaneously. This 
approximation introduces a certain amount of error into the computed H2 abundance. How- 
ever, provided that the timesteps used to evolve the SPH particles in the simulation are 
long compared compared to £ e q,H- and t cqH +, which we have verified is the case for our 
simulations, the size of this error will be negligible compared to that arising from the rate 
coefficient uncertainties that are the subject of this paper. Having made this approximation, 
we are left with four non-equilibrium species, and thus in theory four chemical equations. 
However, since charge conservation implies that 

xr+ + x H + = x e - + x H - , (22) 

while conservation of the total number of hydrogen nuclei implies that 

xr + x H + + xr- + 2x H + + 2xr 2 = 1, (23) 

we can reduce the actual number of equations that we must solve to only two. We choose 
to solve for xr 2 and Xr+, obtaining x e - and rr H from Equations 22-23, but other choices are 
of course possible. Our choice here should make no difference to our final results. 

To solve the chemical equations for a given SPH particle, we make use of a technique 
known as operator splitting. We assume that within the current timestep of the SPH particle, 
we can solve for the density evolution of the gas separately from its chemical evolution. The 
density evolution can then be computed using the same prescription as in standard SPH 
(see Springel et al. 2001 for details), and the updated gas density is then available for use in 
the chemical equations. These are then solved implicitly using DVODE, a freely available 
and well tested double precision ordinary differential equation solver (Brown, Byrne, & 
Hindmarsh 1989). Operator splitting introduces a certain amount of error, as in practice the 
density should vary during the chemical timestep. However, the SPH algorithm naturally 
limits the extent to which the density can change during a single SPH particle timestep, by 
making particles take shorter timesteps in rapidly evolving regions. We therefore expect the 
error introduced by this technique to be small as is discussed further in Section 4. 

In common with other authors, we use a simplified reaction network that does not 
include the chemistry of minor coolants such as HD or LiH. We also neglect any effects 
due to helium chemistry. Neglect of the minor coolants is justified by the fact that H 2 
dominates the cooling of the gas for the all of the temperatures and densities found in our 
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simulations (Flower & Pineau des Forets 2001). At the worst, we may overestimate the 
final temperature of the gas slightly. Neglect of the helium chemistry is also easily justified, 
provided that we assume that the bulk of the helium in the gas is in neutral form, as 
in this case the only reactions involving helium that play any role in determining the H 2 
abundance - the collisional dissociation of H~ and H2 by He - are far less effective than the 
corresponding reactions with H (Abel et al. 1997). So the error in the final H2 abundance 
will be small. If significant amounts of ionized helium are present, then our assumption 
introduces a larger error, since we will underestimate the actual electron abundance in the 
gas, and hence the H~ formation rate. However, even in this case, we would expect the error 
in the final H 2 abundance to be relatively small, owing to the small abundance of helium 
relative to hydrogen. The chemical reactions included in our network are summarized in 
Table 1. In most cases, we also list the source used for the adopted rate coefficient. The 
exceptions are the associative detachment and mutual neutralization reactions, which we 
discuss further below. 

Table 1 lists three photochemical reactions: the photodetachment of H~ and the pho- 
todissociation of H<j~ and H2. To compute the rates for these reactions, we assumed that 
our simulated protogalaxies were illuminated by an external background radiation spectrum 
with the shape of a 10 5 K black-body, as should be typical of the brightest population III 
stars (Cojazzi et al. 2000). We cut off this spectrum at energies greater than 13.6 eV to 
account for absorption by neutral hydrogen in the protogalactic gas and in the intergalactic 
medium. The strength of the background is specified in terms of the flux at the Lyman limit, 
J(v a ) = 10~ 21 J 21 erg s" 1 cm" 2 Hz" 1 sr" 1 . 

If sufficient H 2 forms within the protogalaxy, it will begin to self-shield, reducing the 
effective photodissociation rate. An exact treatment of the effects of self-shielding is compu- 
tationally unfeasible, as it would require us to solve for the full spatial, angular and frequency 
dependence of the radiation field at every timestep. Instead, we have chosen to incorporate it 
in an approximate manner. We assume that the dominant contribution to the self-shielding 
at a given point in the protogalaxy comes from gas close to that point, and so only include 
the contribution to the self-shielding that comes from the nearby H 2 . To implement this 
approximation numerically, we make use of the fact that Gadget already defines a suitable 
local length scale: the SPH smoothing length h, which characterizes the scale over which 
the flow variables are averaged. In Gadget, as in all modern SPH codes, h is allowed to 
vary from point to point within the flow and is automatically adjusted in order to keep the 
mass enclosed within a sphere of radius h roughly constant. Further details can be found 
in Springel et al. (2001). In our calculation of the H 2 column density used to compute the 
degree of self-shielding we include only H 2 which lies within a single smoothing length of the 
point of interest. We justify this approximation by noting that in our simulations, widely 
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separated SPH particles typically move with a significant velocity relative to one another. 
Consequently, the contribution to the total absorption arising from one particle is Doppler 
shifted when viewed from the rest frame of the other particle. If this Doppler shift is large 
compared to the line widths of the Lyman- Werner band transitions, then the effect is to dra- 
matically reduce the extent to which the absorption contributes to the total self-shielding. 
On the other hand, gas close to the point of interest will typically have a much smaller rel- 
ative velocity, and so will contribute far more effectively. Our approximation considers only 
the latter contribution, and assumes that the former contribution is completely negligible. 
In practice, of course, the distant gas is likely to contribute to some non-negligible extent, 
and so we will tend to underestimate the true amount of self-shielding, unless the gas infall 
is highly supersonic. Nevertheless, we believe that our approximation remains useful as it is 
computationally efficient, and also represents an improvement over previous optically thin 
treatments (e.g., Machacek, Bryan, & Abel 2001; Ricotti, Gnedin, & Shull 2002). 

Finally, we assume that ionization from X-rays or cosmic rays is negligible, although 
previous work suggests that even if a low level of ionization is present, it will not have a 
major effect on the outcome of the collapse (Glover & Brand 2003; Machacek et al. 2003). 

A second major modification that we have made to the Gadget code is the inclusion of 
a treatment of radiative heating and cooling. Cooling in our model comes from three main 
sources: electron impact excitation of atomic hydrogen (a.k.a. Lyman-a cooling), which is 
effective only above about 8000 K, rotational and vibrational excitation of H 2 , and Compton 
cooling. Rates for Lyman-a cooling and Compton cooling were taken from Cen (1992), 
while for H2 rovibrational cooling we used a cooling function from Le Bourlot, Pineau des 
Forets, & Flower (1999). In models where an ultraviolet background is present, we include 
the effects of heating from the photodissociation of H 2 , assuming that 0.4 eV of energy per 
photodissociation is deposited as heat (Black & Dalgarno 1977). We also include heating due 
to the ultraviolet pumping of H 2 , following Burton, Hollenbach, & Tielens (1990), although 
this is only important in high density gas (n > 10 4 cm -3 ). 

To incorporate the radiative heating and cooling terms into Gadget, we again use an 
operator splitting technique. In this case, we assume that the change in the internal energy of 
the gas due to pressure work can be computed separately from that due to radiative heating 
and cooling. The former can then be calculated in the same fashion as in the standard 
Gadget code, while the latter can be computed by solving: 

£ = r-A, (24) 

where e is our initial estimate of the internal energy density of the gas, which already includes 
the effects of the pressure work term, T is the radiative heating rate per unit volume and A is 
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the radiative cooling rate per unit volume. We solve this equation implicitly using DVODE 
at the same time that we solve the chemical equations. Again, the use of operator splitting 
introduces some error into the thermal evolution of the gas, but, as before, we expect this 
error to be small (see Section 4 for further discussion). 

Finally, to allow us to represent gas that has collapsed beyond the resolution limit of the 
simulation in a numerically robust manner, we have modified the code to allow it to create 
sink particles - massive, non-gaseous particles, designed to represent dense cores, that can 
accrete gas from their surroundings but otherwise interact only via gravity (Bate, Bonnell, 
& Price 1995). The design and implementation of our sink particle algorithm is discussed 
elsewhere (Jappsen et al. 2005a). 



4. Initial conditions 

As we wish to be able to run a large number of simulations of protogalactic collapse, we 
have chosen to limit the computational cost of each simulation by starting from somewhat 
simplified initial conditions. Since we are not particularly interested (in this paper at least) 
in following the assembly history of the dark matter halo in which the protogalaxy sits, or in 
studying the response of the halo to the cooling of the gas, we chose to model the influence 
of the halo by using a fixed background potential. To construct this potential, we assumed 
that the halo was spherically symmetric, with the density profile of a Burkert halo (Burkert 
1995): 

Pdm(r) = Pdm '° . (25) 
1 + (r 2 /r*) 

We took the central density of the halo to be pdm,o — 0.13 M pc~ 3 and also specified the 
total mass of dark matter in the halo, M dm = 10 7 M . The core radius r c is calculated 
internally by Gadget, based on M dm and the initial redshift z. For the gas, we assumed an 
initially uniform distribution, with an initial density p g , taken to be equal to the cosmo logical 
background density. The initial temperature of the gas was also taken to be uniform, with 
a value T g , the choice of which is discussed below. 

The quantity of gas present in our simulations was taken to be a fraction Q h /Q dm of 
the total mass of dark matter, where f2dm = fi m — fib- We took values for the cosmo logical 
parameters from Spergel et al. (2003), and so fib = 0.047 and fi m = 0.29, giving us a 
total gas mass of M g = 0.047/(0.29 - 0.047) x 10 7 M Q = 1.934 x 10 6 M . In most of our 
simulations, we used 32768 SPH particles to represent this gas, giving each SPH particle a 
mass M part ~ 59 M . In order to properly resolve gravitationally bound clumps (or other 
gravitationally bound structures) in SPH simulations, they must be represented by at least 
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twice as many SPH particles as are used in the SPH smoothing kernel (Bate & Burkert 
1997). In our simulations, our smoothing kernel encompasses approximately 40 particles - 
for reasons of computational efficiency, the number is allowed to vary slightly, but never by 
more than 5 particles - and so our minimum mass resolution is M res ~ 80M part ~ 4720 M . 
To verify that our results are insensitive to the value of M res , we have repeated several of 
our simulations using a larger number of SPH particles, 131072, corresponding to a four 
times smaller value of M res . We find only minor differences (of the order of 10% or less) 
between the temperature and density evolution in these runs and in the lower resolution 
runs discussed below. Since the SPH particle timesteps in these higher resolution runs will 
typically be shorter than those in our lower resolution runs, the size of the error introduced 
by our use of operator splitting will also be different. Therefore, the close agreement of the 
results of these runs with their lower resolution counterparts also helps to reassure us that 
the splitting errors in both sets of runs are small. 

To prevent artificial fragmentation or other numerical artifacts from affecting our results, 
it is necessary either to halt the simulation before the local Jeans mass, Mj falls below M res 
in any part of the simulation volume, or to use sink particles to represent regions where 
Mj < M res . We have chosen the latter course, and so create sinks in regions where the gas 
density exceeds 6 M Q pc -3 using the algorithm described in Jappsen et al. (2005a). This 
corresponds to a hydrogen atom number density n cr i t ~ 200 cm -3 . The Jeans mass, which 
is defined here as 2 

M 3 = ^G-^p g ^cl (26) 

where p g is the gas density and c s is the adiabatic sound speed, has a value at this number 
density of 

(rp \ 3/2 
Took) Mq " (27) 

This is greater than M rcs provided that the temperature of the dense gas exceeds 100 K. As 
we shall see later, the minimum gas temperature reached by dense, gravitationally collapsed 
gas in our simulations is typically no smaller than 150 K and so our simulations remain 
well-resolved up to the point at which sink particles are created. 

Once created, sink particles can accrete gas from their surroundings, provided that 
the SPH particles representing the gas come within a specified radius of the sink particle 
(the accretion radius r acc ), are gravitationally bound to the sink, and satisfy certain other 



2 Note that in this definition, we assume that the self-gravity of the gas dominates on small scales within 
the halo, since the gas is dissipative and so can collapse to much higher densities than the non-dissipative 
dark matter. 
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conditions (see Jappsen et al. 2005a for details). In our simulations, we set r acc = 20 pc. This 
value is chosen to be somewhat larger than Jeans length of the gas at n crit , but is admittedly 
somewhat arbitrary. However, since SPH particles within r acc will only be accreted if they 
satisfy all of the necessary criteria, the outcome of our simulations should be insensitive to 
the precise value chosen for r acc . 

We initialized each of our simulations at a redshift z = 20 and allowed them to run 
for 220 Myr; given our adopted cosmological parameters, this interval corresponds to ap- 
proximately 1.25 Hubble times, with the simulations terminating at a redshift z ~ 11.2. 
Protogalaxies which fail to cool and collapse during this interval are unlikely to get the 
chance to do so thereafter, as the typical interval between major mergers of dark matter 
halos is of the order of a Hubble time (Lacey & Cole 1993). 

Finally, it should be noted that the initial conditions described here are undoubtedly 
highly simplified and are unlikely to properly represent the full hydrodynamics of proto- 
galactic formation. However, they should be more than sufficient for our purposes, as we 
are primarily interested in the difference in the outcomes of simulations run using different 
values for the mutual neutralization and associative detachment rate coefficients, and in this 
respect we would expect the behaviour of these simple models to be a reliable guide to the 
behaviour of more complex, but ultimately more realistic models. 

5. Results 

To investigate how sensitive the cooling and collapse of protogalactic gas are to the choice 
of mutual neutralization and/or associative detachment rate coefficients, we have performed 
five sets of simulations, with initial conditions as summarized in the previous section and 
in Table 2. For each set of parameters listed in Table 2, we performed nine separate runs 
in which the mutual neutralization and associative detachment rate coefficients were varied 
individually. For runs 1-3, we adopted the lower value for the mutual neutralization rate 
coefficient from Dalgarno & Lepp (1987), while for the associative detachment rate coefficient 
we used values of 0.65 x 10~ 9 , 1.30 x 10~ 9 and 5.00 x 10~ 9 cm 3 s _1 respectively, where the 
central value corresponds to the value measured by Schmeltekopf et al. (1967) and the other 
values represent lower and upper bounds on plausible values, as previously discussed in 
Section 2. For runs 4-6, we adopted the mutual neutralization rate coefficient of Croft, 
Dickinson, & Gadea (1999) and varied the associative detachment data as before. Finally, 
for runs 7-9 we repeated this procedure using the mutual neutralization rate coefficient 
of Peterson et al. (1971). For ease of reference, the nine different combinations of rate 
coefficients used are listed in Table 3. 
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Fig. 3. — The density and temperature distribution at the end of run Al. Each point in the 
plot represents the gas density and temperature associated with an individual SPH particle. 

5.1. Cold initial conditions 

Our first set of simulations (runs A1-A9) examined collapse beginning from cold initial 
conditions. In these runs we took initial values for the fractional ionization and H 2 abundance 
from the IGM chemistry model of Standi et al. (1998), and adopted an initial gas temperature 
T g = 12 K, corresponding to the value of the IGM temperature at z = 20 in the absence of 
any form of preheating. 

The state of the gas at the end of run Al is illustrated by the plot of temperature versus 
number density shown in Figure 3. Infalling cold gas has been shock-heated to T > 10 4 K 
within the potential well of the dark matter halo, and some cooling of the highest density 
gas is evident, although this cooling is not yet far advanced. By the end of the simulation, 
the density of the central gas has reached n c ~ 6.5 x 1CT 2 cm -3 and its temperature is 
T c ~ 6900 K (although note that because of the nature of the SPH algorithm, these values 
actually represent a weighted average over gas within a single smoothing length of the center 
of the halo). Its fractional ionization remains fairly low, x e - jC ~ 2.0 x 10 -4 , as does its 
H 2 abundance, x H2 ,c = 1-8 x 10~ 5 , though the latter has increased by nearly an order of 
magnitude over its initial value. The cooling time of this gas is long - approximately 570 Myr 
- although it will decrease significantly once more H 2 forms. 
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The other runs in this set of simulations (runs A2-A9) give very similar results. In 
Table 4, we list the central density, temperature, fractional ionization and H 2 abundance at 
the end of each of these runs. From the table it is clear that while changes in the mutual 
neutralization rate coefficient have only a very minor effect on the final state of the gas, 
changes in the associative detachment rate coefficient have a far more noticeable effect - in 
runs with a high value of the associative detachment rate coefficient, more H 2 is produced, 
and the extra cooling that this provides leads to a lower central temperature and consequently 
a higher central density at the end of the run. For variations of the associative detachment 
rate coefficient within the range that we consider to be plausible, we find a variation of 
approximately a factor of three in the final H 2 abundance. However, the variation in the 
final temperature and density is considerably smaller, approximately 10%, since in none of 
the models has enough H 2 formed to provide particularly effective cooling. 

The fact that the results of these runs are insensitive to the value of the mutual neutral- 
ization rate coefficient is easy to understand if one considers that at a fractional ionization of 
order 10~ 4 , the net rate of destruction of H by mutual neutralization is two to four orders 
of magnitude smaller than the rate of destruction by associative detachment, meaning that 
the latter always dominates. The sensitivity of the results to the choice of associative de- 
tachment rate coefficient comes about in this case due to the influence of a third process, the 
collisional dissociation of H by atomic hydrogen (reaction 15), which has a rate coefficient 
which is competitive with that of associative detachment at high temperature. However, by 
the end of the simulations, the collisional dissociation rate in the central gas has already 
become smaller than the associative detachment rate, and so we would expect its influence 
on future H 2 formation to be slight. 

Ultimately, while the differences we find between the results of the various runs are 
interesting, they do not appear to be particularly important, since in every case the outcome 
of the simulation is essentially the same - the protogalactic gas fails to cool within a Hubble 
time, and so is unlikely to have sufficient time to cool or collapse before the protogalaxy 
merges with a larger structure. 

The fact that gas in these runs fails to cool effectively is not entirely unexpected, given 
our choice of halo mass and initial redshift. Previous work on protogalaxy formation in the 
cold IGM by Tegmark et al. (1997) and others suggests that only gas in halos more massive 
than some critical mass M cr it will cool effectively (reviewed recently by Bromm & Larson 
2004, Ciardi & Ferrara 2005, and Glover 2005). For gas collapsing at z < 20, Tegmark 
et al. (1997) find M crit > 10 7 M , suggesting that cooling in our 10 7 M halos will only 
be marginally effective. We therefore performed two additional sets of simulations, using 
similar initial conditions to runs A1-A9, but starting at z = 30 and z = 40. Appropriate 
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adjustments were of course made to the initial density and temperature of the cold gas. We 
find more evidence for cooling in these runs, in line with the theoretical expectations, but 
variations in the associative detachment and mutual neutralization rate coefficients continue 
to have no more than a small effect, with the largest variation in the central gas density 
being of order 5% and the largest variation in the central H2 abundance being of order 15%. 
We therefore conclude that our z = 20 runs are giving us an accurate picture of the influence 
of the uncertainties in the rate coefficients in this particular scenario. 



5.2. Hot initial conditions 

Since the uncertainties in the rate coefficients do not appear to have a large impact 
on the outcome of protogalactic collapse that begins from cold initial conditions, we chose 
in our next set of runs (B1-B9) to examine an alternative situation, in which we might 
expect the uncertainties to have a greater effect. We took as initial conditions gas which 
was hot (T g = 10 4 K) and fully ionized (x e - = 1.0, x H2 = 0.0). The physical situation that 
these initial conditions are intended to represent is that of a protogalaxy forming within 
what Oh & Haiman (2003) term a 'fossil' Hll region - an Hn region surrounding an ionizing 
source which has just switched off, but the surrounding gas has not yet had time to cool and 
recombine. Prior to cosmological reionization, such regions should be relatively common, 
since the characteristic lifetimes of the likely ionizing sources - massive population 111 stars 
and/or active galactic nuclei - are significantly shorter than the Hubble time. 

In Figure 4, we show the state of the gas at the end of run Bl. The asterisk on the 
right hand side of the plot represents a sink particle of mass M sin k ~ 1.2 x 10 5 M . Since we 
have no information on the true temperature and density distribution of the gas represented 
by the sink particle, beyond the knowledge that its density is greater than 200 cm" 3 , we 
cannot assigned the material within the sink particle to its proper location in this plot. 
We have therefore assigned it a nominal density of 200 cm -3 and temperature of 200 K, as 
these values are broadly representative of the properties of the gas immediately prior to sink 
particle formation. The contrast between this figure and Figure 3 is striking. Far more high 
density gas is present, especially if one takes the material in the sink particle into account, 
and the temperature of the gas is much lower: the minimum temperature reached by the 
end of the simulation in gas which has not become part of the sink is 347 K, while the lowest 
temperature reached in any of the gas during the course of the run is 172 K. Note that the 
reason that the latter value is lower is that we cease tracking the temperature of gas once it 
becomes part of a sink particle; consequently, the minimum temperature appears to increase 
after the creation of the sink particle. Compared to these values, the minimum temperature 
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Fig. 4. — The density and temperature distribution at the end of run Bl. Each point in the 
plot represents the gas density and temperature associated with an individual SPH particle. 
The asterisk represents the nominal temperature and density of the gas comprising the sole 
sink particle immediately before the formation of the sink. 

reached in run Al, which is 6940 K, is a factor of twenty to fourty times larger. 

The reason for this dramatic difference becomes clear if one compares the growth of the 
mass- weighted mean H2 fraction in runs Al and Bl. This quantity is given by 

EjXH 2 ,i 
•^H2,mean 77- 

JVsph 

where xu 2 ,i is the fractional H 2 abundance of the i-th SPH particle and we sum over all 
A sp h SPH particles. Since all of the SPH particles have the same mass, this gives us a 
mass-weighted value. The evolution of XH 2 ,mean during the two runs is plotted in Figure 5. 
Although run Al has a higher initial H 2 abundance than run Bl, it is clear that H 2 forms 
far more rapidly in the latter run, due to the high initial fractional ionization, and it very 
quickly overtakes run Al. Similar differences are seen if instead of comparing the mean H 2 
fraction, we compare the evolution of the H 2 fraction at the center of the density distribution, 
xu 2 ,c- Note also that the apparent decline in XH 2 , m ean at t > 200 Myr is simply an artifact of 
the creation of the sink particle, since we cease tracking the H 2 content of gas that becomes 
part of a sink particle. The much higher H 2 fraction found in run Bl allows the gas to cool 
efficiently to a few hundred K, and consequently enables it to collapse to high density. It is 
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Fig. 5. — The evolution with time of the mass-weighted mean H2 abundance in run Al 
(represented by the dashed line) and run Bl (represented by the solid line). The apparent 
decline in the latter at t — 200 Myr is due to the formation of a sink particle, as explained 
in the text. 

reasonable to suppose that some fraction of this collapsed gas will go on to form stars. 

Given the role played by the high fractional ionization in boosting H 2 formation in this 
set of runs, we might expect the uncertainties in the associative detachment and mutual 
neutralization rate coefficients to have far more effect here than they did in runs A1-A9. 
Indeed, this is what we find. In Figures 6-8, we show how the central H 2 abundance, central 
gas density, and central gas temperature vary over the course of runs B1-B9, while in Figure 9 
we plot the evolution with time of f cc , the fraction of the gas which is cool and condensed, 
defined here as the fraction of the total gas which has a temperature T < 500 K and a density 
n > 100 cm~ 3 . 

Figure 6 demonstrates that the differences in the associative detachment and mutual 
neutralization rate coefficients used in these runs lead to large differences in the resulting 
H 2 abundances. For instance, the difference between the value of ih 2)C in runs B3 and B7 
at t — 150 Myr is greater than an order of magnitude; and while there is some indication 
that the value of xh 2 ,c in the different runs begins to converge at late times, the differences 
between the runs remain significant at the end of the simulations. It should also be noted 
that the sense of the differences is precisely what we would expect, based on our previous 
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Fig. 6. — The evolution with time of the central H 2 abundance found in the gas in runs Bl- 
B9. From bottom to top, the dotted lines correspond to runs B7, B8, and B9 respectively, 
the dashed lines to runs B4, B5, and B6 and the solid lines to runs Bl, B2, and B3. For 
clarity, we only plot the evolution up until the point at which a sink particle forms (or until 
the end of the run, if no sink forms). 
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Same as Figure 6, but for the central density of the gas. 
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Fig. 8. — Same as Figure 6, but for the central temperature of the gas. Note that from 
bottom to top, the solid lines now correspond to runs B3, B2, and Bl respectively, the 
dashed lines to runs B6, B5, and B4 and the dotted lines to runs B9, B8, and B7. 
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Fig. 9. — The evolution with time of / cc , the fraction of cool, condensed gas. The labeling 
of the lines is as in Figure 6. Note that in runs B4, B7, and B8, f cc remains zero until the 
end of the run. 
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discussion of the chemistry: H 2 formation is more efficient in runs with a large value of the 
associative detachment rate coefficient than in runs with a smaller value, but is less effective 
when the mutual neutralization rate coefficient is large than when it is small. 

Figures 7 and 8 follow the evolution of the density and temperature at the center of the 
protogalaxy and demonstrate the effect that these differences in H 2 abundance have on the 
thermal and dynamical evolution of the gas. For the first 70Myr there is very little difference 
between the runs. In each case, the gas initially cools to T ~ 4000K through a combination of 
Lyman-a cooling (at T > 8000 K) and Compton cooling, before gradually being re-heated by 
adiabatic compression in the gravitational potential well of the dark matter. At later times, 
the gas begins to cool more strongly through H2 ro-vibrational emission, and the ensuing 
loss of pressure support leads to the gas rapidly collapsing to high densities. However, the 
time of onset of this phase of cooling and collapse is sensitive to the H 2 abundance of the gas, 
and so occurs significantly later in runs with small H 2 abundances than in runs with high 
H 2 abundances. The physical reason for this dependence is clear: since H 2 is the dominant 
coolant in the gas at t = 70 Myr, the cooling time of the gas scales roughly as t coo i oc x^, 
and so protogalaxies in which the gas forms less H 2 naturally take longer to cool. 

One important consequence of this is that the runs display a wide variation in the 
amount of cold dense gas available for star formation. This is illustrated by our plot of f cc 
in Figure 9, which shows that the fraction of cold, dense gas varies from 0% (in runs B4, B8 
& B9) to almost 25% (in run B3). Given a sufficiently long time to evolve, it is likely that 
a significant amount of cold gas will eventually accumulate in the runs with less efficient 
H 2 formation. However, we must recall that unlike our idealized model protogalaxies, real 
protogalaxies are not isolated systems, and do not have an unlimited time in which to evolve 
before being disrupted by mergers or other external events (such as nearby galactic outflows; 
see Thacker, Scannapieco, & Davis 2002). Therefore, our determination of whether or not 
a given halo can cool quickly enough to form stars may depend on our choice of associative 
detachment and mutual neutralization rate coefficients. 



5.3. Cooling and collapse with a UV background 

So far, we have assumed that the effects of any external UV background are negligible. 
In practice, this is unlikely to be the case - by z = 20 a considerable UV background may 
already have developed (Haiman, Abel, & Rees 2000; Glover & Brand 2003). Indeed, if 
cosmological reionization is to occur somewhere in the redshift range z re i on = 17 ± 5, as is 
suggested by the WMAP polarization results (Kogut et al. 2003), there must already be a 
fairly strong background in place. To explore how the presence of a UV background may 
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influence our conclusions, we have run several sets of simulations using the same initial 
conditions as runs B1-B9, but in which the strength of the UV background was varied as 
listed in Table 2. 

In Figures 10-12, we show how the central H 2 abundance varies with time in runs Cl- 
C9, D1-D9 and E1-E9, which correspond to UV field strengths specified by J 2 i = 1CT 2 , 
3 x 10~ 3 and 10~ 3 respectively. From Figure 10, it is apparent that the presence of a strong 
UV background significantly enlarges the difference in outcomes between runs with different 
associative detachment and mutual neutralization rate coefficients. For instance, the central 
H2 abundance at the end of run C3 differs from that at the end of run C7 by almost a 
factor of 100, whereas runs B3 and B7 differ in this respect by no more than a factor of 4. 
However, from Figures 11 & 12 we see that as the field weakens, the differences become less 
pronounced; indeed, the results of runs E1-E9 are close to those of runs B1-B9. 

The reason for the large spread in outcomes that we see in runs C1-C9 becomes clearer 
once we consider that photodissociation by the imposed UV background will set an upper 
limit on the H 2 abundance in optically thin gas, given by the equilibrium value at which 
formation balances photodissociation: 

R 

^H 2 ,eq ^ "5^H. (29) 

-Kpd 

Here R p d is the photodissociation rate per H 2 molecule, and Ru 2 is the rate of H 2 formation 
per H atom. To compute the latter, we note that the dominant contribution comes from H 2 
formation via the H~ pathway, and so we can approximate i? H2 by the product of the H 
formation rate per H atom and the fraction of H~ that successfully forms H 2 (rather than 
being destroyed by mutual neutralization). Therefore 

r> 1 ^ad^H /OA \ 

Rh 2 = k H ~n e - x — , (30) 

^ad^H + k mn n H + 

where ku- is the rate coefficient for the formation of H~ by radiative association, and where 
A; ad and /c mn are the associative detachment and mutual neutralization rate coefficients re- 
spectively. In gas with a high fractional ionization, /c mn %+ & a d^H, and so this expression 
reduces to 

Ru 2 ^ ^H-T^^^H, (31) 

fc mn n H + 

which can be further simplified to 

^ad 

Rn 2 - ^h-t— ^h, (32) 

provided that n e - ~ n H +. After 100 Myr of evolution, the gas at the center of our simulated 
protogalaxies in runs C1-C9 has a temperature of approximately T = 7100 K and an atomic 
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Fig. 10. — The evolution with time of the central H2 abundance in runs C1-C9. The ordering 
of the lines is the same as in Figure 6. 
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Fig. 11. — Same as Figure 10, but for runs D1-D9. 
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Fig. 12. — Same as Figure 10, but for runs E1-E9. 

hydrogen number density of approximately % = 2.6 x 1CT 2 cm -3 ; the variation in these 
values between the different runs is about 5-10%. Using these values for % and T, we find 
that Aft- = 3xl0- 15 cm 3 s _1 (Wishart 1979), and hence that Ru 2 = 3 x 10 _15 (A; ad /A; mn )nHS- 1 . 
Now, since R pd = 1.3 x 10~ 14 s _1 , given our assumed UV field strength of J 2 i = 10~ 2 (Draine 
& Bertoldi 1996), it is easy to show that 



where the size of the term in brackets varies from 5 x 10~ 3 to 6 x 10 _1 depending on the 
values chosen for the rate coefficients. The equilibrium H 2 abundance in optically thin gas 
therefore varies from 3.1 x 10~ 5 to 3.6 x 10~ 3 , depending on our choice of chemical rate 
coefficients. If we take the effects of self-shielding into account, we will obtain slightly larger 
values, but if anything the spread in outcomes will be even greater, since gas with a large 
H 2 abundance can self-shield far more effectively than gas with a small H 2 abundance. 

If we now compute the cooling timescale for this gas, which is given by 



where n tot is the total number density of particles (including electrons), and where A is the 
cooling rate, computed using the SPH code in units of erg cm -3 s _1 for a range of different 
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H 2 abundances (see Figure 13), then we can determine how much H 2 is required in order to 
cool the gas before the end of the simulation. From Figure 13, we see that an H 2 abundance 
Xh 2 > 10~ 4 is required (as here we are at 100 Myr and the simulations run for 220 Myr). 
Therefore, for those of our runs in which XH 2 ,eq is larger than this value, we would expect 
the gas to cool efficiently, while in runs in which XH 2 ,eq is not sufficiently large, we would 
expect the gas to cool slowly, if at all. If we examine the evolution of the temperature 
of the gas at the center of the protogalaxies simulated in these runs, which we plot in 
Figure 14, we find that our results are consistent with these expectations: in runs C6, C8 
& C9, where XH 2 , e q < 10~ 4 , very little cooling occurs, while in the remainder of the runs, 
which have iH 2 ,cq > 10~ 4 , far more cooling occurs. It is also clear from the plot that the 
final temperature reached in the runs that cool depends strongly on the value of iH 2 ,eq- The 
differences between the thermal evolution of the gas in the various runs are also reflected 
in the dynamical evolution, as can be seen from Figure 15. This example therefore gives a 
particularly striking demonstration of the need for better rate coefficient data, since without 
it we cannot determine which of the various outcomes illustrated in Figures 10, 14 & 15 is 
actually the correct one. 

It is also interesting to compare the length of the cooling time in these runs with the 
free-fall time, given by 



where p tot is the total matter density (i.e., the sum of the gas density and the dark matter 
density). At this point of 100 Myr in the simulation, the main contribution to p tot comes 
from the dark matter. Evaluating tg at the center of the dark matter halo, we find that 
is has a value of approximately 20 Myr. From Figure 13, it is clear that for the gas to 
have t coo \ < tg, it must have an H 2 abundance x^ 2 > 10~ 3 . Since this value is an order of 
magnitude larger than the abundance required to allow cooling by the end of the simulation, 
i.e., within what is essentially a Hubble time tn, it is possible for there to be protogalaxies 
in which tg < t coo \ < tn- In these protogalaxies, the timescale on which the gas collapses will 
be determined by t coo i, rather than tg, and so will be directly sensitive to the amount of H 2 
within the gas. Therefore, predictions that we make concerning the dynamical evolution of 
such protogalaxies will be highly sensitive to our choices for the mutual neutralization and 
associative detachment rate coefficients, as can already been seen from Figure 15. 

If the strength of the UV background is reduced, as in runs D1-D9 or E1-E9, then the 
value of a;H 2 ,eq increases, and so the uncertainty in the rate coefficients has less influence, 
as it is easier for the gas to form enough H 2 to provide efficient cooling. Nevertheless, 
the differences between the runs remain significant, as can be seen from the plots of H 2 
abundance versus time in Figures 11 and 12, or the plots of temperature versus time shown 
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H„ abundance 



Fig. 13. — The cooling time of gas as a function of its H2 abundance, computed assuming a 
temperature T = 7100 K and density % = 2.6 x 10~ 2 cm -3 . These values are representative 
of the central temperatures and densities found within the halos simulated in runs C1-C9. 
Values in the individual runs differ slightly from these values, but never by more than 5-10%. 

in Figures 16 and 17. 

One final comparison that we can make between the runs is to examine how sensitive 
f cc is to the values of the associative detachment and mutual neutralization rate coefficients, 
and to the strength of the UV background. In Table 5 we list the value of f cc at the end 
of the simulation for all of the runs in sets B, C, D and E. The trend of increasing f cc 
with decreasing field strength is clear, and is easily understood given our discussion above. 
However, the table also demonstrates the degree of variation in f cc that is introduced by the 
uncertainties in the rate coefficients. In particular, we see that the lack of reliability in the 
atomic data introduce an uncertainty of almost an order of magnitude into the value of the 
UV field strength that is required to prevent gas cooling and collapse. 

6. Discussion 

The results that we have presented in this paper serve as a demonstration of the potential 
impact that the existing uncertainties in the values of the associative detachment and mutual 
neutralization rate coefficients will have on our ability to make accurate predictions regarding 
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protogalactic cooling and star formation. 

For protogalaxies forming from cold gas with a low initial fractional ionization, the 
impact is small - some uncertainty is introduced into the predicted H 2 abundance, but not 
enough to significantly alter the dynamical evolution of the gas. If we were to continue our 
simulations for a longer period, then we would expect H 2 formation in the cooling, dense gas 
at the centre of the protogalaxy to proceed at a very similar rate in all of our runs. 

For protogalaxies forming from hot, highly ionized gas, however, the impact is substan- 
tial. The uncertainties in the chemical rate coefficients are responsible for producing a large 
uncertainty in the predicted H 2 abundance, which in turn introduces a large uncertainty 
into the cooling rate of the gas. Since the timescale on which the gas collapses is set by the 
larger of t coo i and tg, and since, in most protogalaxies, t coo i is initially longer than tg, any 
uncertainty in the cooling rate has a direct impact on the dynamical evolution of the gas. 
Indeed, in the examples studied in this paper, our choice of rate coefficients has a great influ- 
ence on the conclusion that we draw regarding whether or not the gas can cool and collapse 
within the lifetime of the protogalaxy (which will be of the order of a Hubble time). The 
presence of an ultraviolet background only serves to exacerbate this problem, as it increases 
the sensitivity of the predicted H 2 abundance to our choice of chemical rate coefficients. 

Ultimately, the only way to remove these uncertainties will be to obtain more accurate 
rate coefficients for the relevant associative detachment and mutual neutralization processes. 
Most important for this will be new laboratory measurements for each reaction at cosmo- 
logically relevant collision energies supplemented by further theoretical calculations. 
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Table 1. A list of all the reactions included in our model of primordial gas chemistry. 

Reaction Number Reaction Reference 

1 H + e-^H- + 7 Wishart (1979) 

2 H + H -»• H 2 + e" See text 

3 H + H+ -»■ H+ + 7 Ramaker & Peek (1976) 

4 H + fk| -»• H 2 + H+ Karpas, Anicich, & Huntress (1979) 

5 H- + H+ -> H + H See text 

6 H- + 7 -> H + e" Wishart (1979) 

7 H+ + c- -> H + H Schneider et al. (1994) 

8 H 2 + H+ -»■ H+ + H Savin et al. (2004) 

9 H 2 + e^ -> H + H + c- Stibbe & Tennyson (1999) 

10 H 2 + H -> H + H + H Mac Low & Shull (1986) 

11 H 2 + 7 ^H + H Draine & Bertoldi (1996) 

12 H + e- ^H+ + e-+e" Janev et al. (1987) 

13 H+ + c- -> H + 7 Ferland et al. (1992) 

14 H~ + c~ -> H + e~ + e~ Janev et al. (1987) 

15 H~ + H -»• H + H + e" Janev et al. (1987) 

16 H- + H+ -> H+ + e~ Poulaert et al. (1978) 
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Table 1 — Continued 



Reaction Number Reaction 



Reference 



17 



H+ + 7 -> H + H+ Dunn (1968) 



References. - References are to the primary source 
of data for each reaction. Photochemical reactions as- 
sume an incident spectrum corresponding to a modified, 
diluted black-body, as described in the text. 



Table 2. Initial conditions for the simulations. 



Run T g (K) x e - xh 2 J21 



A 
B 
C 
D 
E 



12 
10 4 
10 4 
10 4 
10 4 



2.2 x 1(T 4 
1.0 
1.0 
1.0 
1.0 



2.4 x 10~ 6 
0.0 
0.0 
0.0 
0.0 



0.0 
0.0 

10- 2 

3 x 10~ 3 
10~ 3 
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Table 3. The values of the mutual neutralization and associative detachment rate 
coefficients used in our runs. We list here all nine possible combinations, together with the 
number used elsewhere in the paper to refer to each combination. 



Set Rate coefficient (cm 3 s x ) 

Mutual Neutralization Associative Detachment 



1 


7 x 10- 


7 T - 


-1/2 


0.65 x 


10" 


-9 


2 


7 x 10 _ 




-1/2 


1.30 x 


10" 


-9 


3 


7 x 10- 


7 T - 


-1/2 


5.00 x 


10" 


-9 


4 


2.4 x 10- 6 (l + 5 


x l^T)T- 1 / 2 


0.65 x 


io- 


-9 


5 


2.4 x 10- 6 (1 + 5 


x 10- 5 T)T- 1 / 2 


1.30 x 


10" 


-9 


6 


2.4 x 10- 6 (1 + 5 


x 10~ 5 T)T- 1 / 2 


5.00 x 


10" 


-9 


7 


5.7 x lO^T- 1 / 2 
-9.2 x lO^T 1 / 2 


+ 
+ 


6.3 x 10~ 8 

4.4 x 10- 13 T 


0.65 x 


io- 


-9 


8 


5.7 x lO^T- 1 / 2 
-9.2 x lO^T 1 / 2 


+ 
+ 


6.3 x 10- 8 

4.4 x 10- 13 T 


1.30 x 


io- 


-9 


9 


5.7 x lO^T- 1 / 2 
-9.2 x lO^T 1 / 2 


+ 
+ 


6.3 x 10- 8 

4.4 x 10- 13 T 


5.00 x 


io- 


-9 



References. - The mutual neutralization rate coefficients are taken 
from Dalgarno & Lepp (1987), Croft, Dickinson, & Gadea (1999) and 
Peterson et al. (1971) respectively. For the associative detachment rate 
coefficients, the central value is taken from Schmeltekopf et al. (1967); 
the other values represent lower and upper bounds on plausible values. 
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Table 4. Physical state of the densest gas at the end of runs A1-A9. 



Run 


rri 1 1 r\'A T/ 

Jc/lO K 


/1 n— 2 —3 

n c /10 cm 


/I A — 4 


/ 1 n— 5 
^H 2 ,c/10 


Al 


6.94 


6.47 


1.98 


1.83 


A2 


6.66 


6.75 


1.98 


3.08 


A3 


6.30 


7.30 


1.98 


5.46 


A4 


7.72 


6.54 


1.98 


0.66 


A5 


7.44 


6.71 


1.98 


1.24 


A6 


6.35 


7.32 


1.98 


5.06 


A7 


6.92 


6.48 


1.96 


1.85 


A8 


6.66 


6.79 


1.96 


3.00 


A9 


6.30 


7.25 


1.97 


5.17 



Table 5. The value of the cold, collapsed gas fraction, / cc , at the end of the simualtions. 



Run 


J21 = io~ 2 


J 21 = 3 x 10" 3 


J21 = 10- 3 


J21 = 0.0 


1 


0.00 


0.00 


0.08 


0.12 


2 


0.00 


0.08 


0.14 


0.17 


3 


2 x 10~ 4 


0.15 


0.19 


0.23 


4 


0.00 


0.00 


0.00 


0.00 


5 


0.00 


0.00 


0.005 


0.08 


6 


0.00 


0.10 


0.15 


0.18 


7 


0.00 


0.00 


0.00 


0.00 


8 


0.00 


0.00 


0.00 


0.00 


9 


0.00 


0.00 


0.07 


0.11 



